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We discuss convergence and coupling of Markov chains, and present general relations between the 
transfer matrices describing these two processes. We then analyze a recently developed local-patch 
algorithm, which computes rigorous upper bound for the coupling time of a Markov chain for non- 
trivial statistical-mechanics models. Using the "coupling from the past" protocol, this allows one 
to exactly sample the underlying equilibrium distribution. For spin glasses in two and three spatial 
dimensions, the local-patch algorithm works at lower temperatures than previous exact-sampling 
methods. We discuss variants of the algorithm which might allow one to reach, in three dimensions, 
the spin-glass transition temperature. The algorithm can be adapted to hard-sphere models. For 

O two-dimensional hard disks, the algorithm allows us to draw exact samples at higher densities than 
previously possible. 
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The Monte Carlo method is a fundamental computational tool in science. Its goal is to sample configurations 
H 1 in a given state space from a probability distribution n(x). This can usually not be achieved directly for multi- 
dimensional distributions. Markov-chain Monte Carlo methods [1, 2] overcome this problem by generating configura- 
tions xq, X\, x^i ■ ■ ■ starting from an initial configuration xq which belongs to a simpler distribution ir° (often a fixed 
' initial condition, or some ad-hoc random choice). Configurations Xt are then generated from Xt~i according to a 
stochastic algorithm which guarantees, as time moves on, that 7r* departs from the initial condition and converges for 
t — > 00 towards the equilibrium distribution 7r°° = tt. The Markov-chain approach can be implemented for arbitrary 
distributions tt, using for example the Metropolis and the heat-bath algorithms. For many applications, enormous 
effort has gone into designing fast algorithms for which one reaches n (x) ~ n for reasonable running times t. In this 
paper, we are concerned with a related problem: rather than to find the fastest algorithm for a given problem, we are 
interested in quantifying the speed of a given Markov-chain algorithm. This is, we want to prove after which time 
£j I t the sample Xt is equilibrated. It then reflects the equilibrium distribution and no longer the initial configurations. 
1 — In many practical applications, it is extremely difficult to decide from within the simulation whether it has indeed 
equilibrated [2-4]. Instead, one must validate the simulation results with other approaches, from exact solutions to 
experimental data. The correct characterization of the convergence towards equilibrium from within the simulation 
■ has remained a serious conceptual and practical problem of the Monte Carlo method. 

\ From a fundamental viewpoint, the problem of rigorously proving convergence of a simulation was solved, at least 
in principle, through a paradigm called "exact sampling" , which allows to generate, with Markov chains, samples x 
directly from the equilibrium distribution n without any influence of the initial configuration [5]. In practice, however, 

. it has not been possible to implement exact sampling for many complicated problems, as for example disordered 
t-H ' systems, for which standard methods of evaluating equilibration times fail. 

0^ . The reason for this difficulty is as follows: Exact sampling proves for a given Markov-chain simulation that the 
correlation of the initial configuration with the configuration at time t strictly vanishes. This is done by showing 
explicitly that all possible initial configurations x to yield the same output under coupled Monte Carlo dynamics. In 
many simple models, one can prove this coupling property indirectly. In general, however, one must indeed survey 
the entire configuration space. This is usually too complicated to be achieved. 

We recently developed a local-patch algorithm[6] which indeed monitors the entire configuration space of complicated 
systems, even for very large sizes. The approach uses local information, concentrated on so-called "patches" . The scale 
of these patches increases during the simulation. Information on patches can then be combined for the entire system 
to provide a crucial upper bound for the (global) coupling time, and to generate an exact sample. The algorithm 
was demonstrated to work for spin glasses at lower temperatures than previous methods [7, 8], even though the 
physically interesting regime has still not been reached yet. The local-patch algorithm is quite general: in addition 
to spin glasses, we implement it in this paper for hard disks and improve on previous results [9, 10]. The successful 
application of exact sampling to hard-sphere systems is remarquable because the configuration space is continuous so 
that, naively, its complete survey appears out of reach. 
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A. Transfer matrix 

A Markov chain is fully characterized by the so-called "transfer matrix" of transition probabilities between any 
two configurations k and I. As will be illustrated shortly (Section II) in a specific example, the largest eigenvalue of 
the transfer matrix is Xi — 1 and the corresponding eigenvector ^ describes the equilibrium state. The convergence 
towards equilibrium is governed by the spectrum of the transfer matrix and by the overlap of the eigenvectors 
with the initial configuration: 

Ax)=ir(x)+ J2 <*fck°> (1) 

*fc,A fc <l 

In the limit of infinite simulation time, the second-largest eigenvalue determines the exponential convergence of the 
probability distribution towards equilibrium. This eigenvalue sets a time scale 

Tcorr = 1/| log (A 2 ) |, 

and the convergence is as 

7T* — 7T CX exp ( — i/ T corr) for t — > OO. (2) 

The rigorous determination of convergence properties of Markov chains has been undertaken in many cases, from urn 
models to card-shuffling (see [11]), diffusion processes, and many more (see [12]). Efficient algorithms, as for example 
the bunching method [2] are commonly used to perform an empirical error analysis of Monte Carlo data in more 
complicated cases, where rigorous calculations are out of the question. However, these methods are not failsafe. In 
practice, it is often difficult to extract r corr from the large number of physically relevant time scales. In disordered 
systems, for example, there is often no reliable way to ascertain that the simulation has run long enough, and r corr 
may be much larger than assumed (see e.g. [2], Sect. 1.5). 

B. Loss of correlation and exact sampling 

In the limit of infinite times t — ► oo, the Markov chain converges towards the equilibrium distribution, and the 
positions Xt become independent of the initial condition. The loss of correlation with the initial condition is evident 
for Markov chains that couple, that is, which for each possible initial condition xq produce the same output x t . In 
many cases of interest this happens after a finite global "coupling time" t > i C ou P , which depends on the realization 
of the Markov chain. Propp and Wilson [5] realized that this coupling property allows one to draw "exact" samples 
from the distribution n. 

This approach, called "coupling from the past", eliminates the problem of analyzing the convergence properties. 
However, to establish that a Markov chain has coupled, the entire state space of the system must be supervised. This 
was believed infeasible except for special problems where the dynamics conserves a certain (partial) ordering relation 
on configurations. A partial order is conserved in heat-bath dynamics of the ferromagnetic Ising model, whereas the 
frustration in the spin-glass model foils this simplification. 

II. COUPLING AND CONVERGENCE IN A ONE-DIMENSIONAL MODEL 

We first discuss convergence and coupling for a Markov chain describing the hopping of a single particle on a simple 
iV-site lattice with periodic boundary conditions (see Fig. 1). In one time step, the particle hops with probability | 
from one site to its two neighbors: 

Vk^k+i =Pk^k-i = 1/3 (if possible). (3) 

In addition, we have n — Pjv->i = 1/3. 

The equal hopping probabilities imply via the detailed balance condition 

KkPk^i = npi^k (4) 

that the stationary probability distribution -k^ — 1/N of this problem is independent of k. 
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FIG. 1: A Markov chain on a five-site lattice with periodic boundary conditions. The particle hops from a site k towards its 
neighbors with probability 1/3 each. 



This system's Monte Carlo algorithm is encoded in the N x N transfer matrix T 11 : 



T 1 ' 1 = {p(i^j)} = - 
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The eigenvalues of T 1 - 1 are A*' 1 = § (l + 2 cos 2( V )7r ) ,& = !,•■ ■ ,Int[iV/2] + l (with multiplicities) that is, for N = 5, 



{1, ^p}. The largest eigenvalue, Aj 



i,i 



1 corresponds to the conservation of probabilities. By construction, 

1,1 



it is associated with the equilibrium solution (m, . . . , 7Tjv) — • • • , j?)- The second-largest eigenvalue is A 2 ' . For 
,1,1 _ l+y^ 

which vanish as [A 2 ' ]* = cxp [— t/r corr ], with 



N = 5 we have A2' 1 = 1+ 6 V& = 0.539. This eigenvalue controls the long-time corrections to the stationary solution, 



1/1 log (A^)|. 



We note that the time scale r corr only describes the asymptotic behavior of the correlation. The calculation of the 
time t at which the probability distribution 7T* itself is within a suitably chosen e of a the equilibrium distribution 7r 
is more involved (see, for example, [11, 12]). 



A. Coupling 

As illustrated in Fig. 2, the Monte Carlo algorithm can be formulated in terms of random maps. In our example, 
this means that instead of prescribing one move per time step, as in Fig. 1, we now sample moves for all times t 
and all sites k, in such a way that the dynamics of a single particle again satisfies the detailed balance condition of 
Eq. (4). The most natural implementation of this approach is illustrated in Fig. 2: arrows are chosen independently 
for all times t and all sites k. At time to, for example, the particle should move down from sites 1, 3, 4 and 5 and 
straight from site 2. We can now check the outcome of the Monte Carlo calculation. In the example of Fig. 2, from 
time to + 10 on, all initial configurations of the single particle yield the same output. This is remarkable because, 
evidently, at this time the initial conditions are completely forgotten. 

The coupling time i C oup is a random variable (t C oup = 10 in Fig. 2) which depends on the realization of the full 
Monte Carlo simulation from time to onwards (until coupling has been reached). The independence of random maps 
on different time steps implies that the probability for not coupling vanishes at least exponentially fast in the limit 
t — > 00. 

Under the random-map dynamics, an initial state with N particles eventually evolves into a state with one particle 
(in later sections, spin-glass configurations will take the place of the single-particle positions). More generally, a state 
with k configurations can evolve at each time step into a state with k' < k configurations. Figure 2 displays a sequence 
of random maps and illustrates the associated time-forward search of the coupling time. 



4 




FIG. 2: Extended Monte Carlo simulation on N = 5 sites. Trajectories from all possible initial configurations at t — to are 
indicated. They "couple" at t = io + icoup- The coupling time (here t C ou P = 10) depends on the realization of the Markov chain. 



This extended Monte Carlo dynamics on fc-configuration states can again be described by a transfer matrix: 

... \ 



T 3 < 3 



V 







(6) 



where the block T k ' 1 (of sizes (JTj x (ff)) concerns all the processes which lead from a state at time t with k 
configurations to a state with I < k configurations at time t + 1. The upper left block of this matrix, T 1 ' 1 , is the 
original matrix from Eq. (5). As an example, we find from Eq. (3) the following elements of this transfer matrix: 

T fw {| oo.o.) -» I o.o.o}} = 1/9 
T lw {\ oo.o.} -» I ooo.o}} = 1/9 
T fw {| oo.«.) ^ I ooo.o}} = 1/27, 

etc. The matrix T describes a physical system with variable particle number (from 1 to N) and a space comprising 
2 N — 1 states, the number of non-empty states in this new simulation (for a problem of N spins, the number of 
configurations is 2 N and the total number of fc-configuration states (states with k configurations) is 2 2 — 1). 
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FIG. 3: The exact probability that the Markov chains have not coupled by time t (computed by repeated application of the 
forward transfer matrix), compared to the time-scales Tcorr, T C oup , and t* (One-dimensional diffusion model with N = 5). 



The "forward" transfer matrix T allows us to compute the coupling probabilities as a function of time in Fig. 3. 
The matrix T fw is block-triangular in the number of particles (k,l), with the (1, 1) block given by T 1 ' 1 . Therefore, 
all the eigenvalues of T 1,1 are also eigenvalues of T fw . In particular, the largest eigenvalue of T tw is again A^ w = 1, 
with corresponding right eigenvector |(|ioooo) + ---+|oooo»)). The second-largest eigenvalue of T fw belongs to 
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the (2, 2) block and leads to the time scale of the coupling, r COU p- It is given by A| w = 0.838, larger than A MC = A2' 1 . 
This second-largest eigenvalue A| w governs the coupling probability V{t COVLV ) for large times. It follows from the block- 
triangular form of the forward transfer matrix T fw that the time scales satisfy r coup > r corr . A general argument allows 
us to better understand this result: for any running time t we may separate all the Markov chains into those that have 
already coupled and those that have not. Only the non-coupled chains (whose number vanishes as exp (— t/r coup )) 
contribute to connected correlation functions: 

cxp (-Vw) oc (O(t)O(0)) = ]T V(t coup > t, <J t , a )O(a )O{a t ) 

config. a t , a a 

+ J2 K°MO(<7 ) P(*«mp<*,0t)0(<7t). 
config. cto config. at 

Here, O is an observable whose mean value is zero and a t the configuration of the system at time t. For the chains 
which have coupled by the time t, at does not depend on ao, and the contribution to the correlation function vanishes. 
For the other chains, we find 

cxp (-t/Vcorr) = ^2 V(t coup > t, (T t ,a )0((7o)0((T t ) oc exp (-t/r coup ) exp (-t/r*) , (7) 

config. a t , cr 

where we suppose that even the non-coupled chains converge towards equilibrium on a time scale r*, and use that 
the probability for a chain not to have coupled behaves as exp (— i/r coup ) in the long-time limit. Equation (7) shows 
that the difference between r coup and T corr is caused by the convergence taking place within non-coupling chains: 



7"corr 7"coup T 

This relation is illustrated in Fig. 3 for the one-dimensional diffusion model with N — 5. A later figure, Fig. 8, will 
illustrate for the case of spin glasses the split between the general spin-spin correlation function and that same object 
computed for non-coupling chains only. 



B. Forward and backward coupling 

The probability distribution of coupling times in the forward direction can be obtained from the transfer matrix 
T fw as we discussed in Section II A. Here we analyze the distribution of coupling times in the backward direction 
for the application of the "coupling from the past" protocol which, as we will see, is the same as the one in the 
forward direction. The backward coupling process leads to a generalized transfer matrix, T bw , which again describes 
an extended Monte Carlo simulation. 




FIG. 4: Extended simulation on N = 5 sites. The outcome of this simulation, from t = —00 up to t = 0, is k = 2. It can be 
obtained by backtracking from time t = to — t^ up , or by forward simulation from any to < ~t^ up , through the indicated 
trajectories. Backtracking from sites 1,3,4 and 5 leads to dead ends. 

We consider a hypothetical simulation which has run since time t = —00 up to time t = (see Fig. 4). It follows 
from the discussion in Section II A that the simulation has coupled. Furthermore, because of the infinite separation 
between the infinitely remote initial condition and the final one, the resulting configuration (at t = 0) is in equilibrium. 
But it remains to be seen which one of the five configurations at t = is generated. In the example of Fig. 4, a 
one-step backtrack to time t = — 1 allows us to see that the configuration at t = can be neither k = 1 nor k = 3 
nor k = 5. Likewise, for N — 1 output positions on the ^V-site ring this back-propagation leads to a dead end, and 



6 



only a single position yields a full set of possibilities at some time — t%£ up . Thus to find the output configuration of 
the simulation, one is interested in the first time in the past for which the simulation couples, that is one searches 
the "backward" coupling time. The implementation of this backward simulation, as defined for many particles, can 
again be described by a transfer matrix. For any distribution of arrows, any occupied site fc at time t propagates 
its occupation back to all the sites at time t — 1 which have arrows pointing towards fc. The matrix element of T bw 
between two states is given by the statistical weights of all the arrows connecting the two states. For example, we 
find for T bw 



T bw {| oo.oo) - 


* | o • • o o)} 


= 4/81 


T bw {| ••o..} - 


» | • o o • •)} 


= 4/81 


T bw {| oo.o.}} - 

S v ' 


> o o o • o)} 


= 2/27. 


time t 


time t + 1 





This is a non-trivial variant of the forward simulation as, for example, the matrix T bw is not block-triangular as T , 
but it is particle-hole symmetric (as we see in the above example). 

Formally, a random map / (here a set of arrows) associates configurations which are connected under the Monte 
Carlo dynamics. A fc-configuration state \x) is by definition a set of configurations. At time t the state \x) can be 
associated with the state \y) at time t + 1 via the forward matrix if and only if \y) is the (set) image of \x) by an 
allowed mapping / (i.e. \y) = /(|a;))). The same holds for the backward matrix but |a;) must be the reciprocal (set) 
image of \y) (i.e. |a;) = / _1 (|y)))- The backward transfer matrix T bw manifestly differs from the forward matrix T fw . 
However, we construct explicitly in Appendix A the similarity transform that maps T fw onto T bw . This means that 



(see Appendix A). The similarity transformation P associates a fc-configuration state \x) with the sum of states that 
share at least one configuration with \x), included itself. The spectrum of the backward transfer matrix thus agrees 
with the one of the forward matrix and the distribution of coupling times V(t coup ) is identical for backward and 
forward dynamics. This result is natural because the probabilities for not coupling for t time steps are identical in 
both forward and backward direction: V(t^ up > t) = V(t l ™ up > t) (see [9]). The probability distribution V(t = i b ™ up ) 
measures the weight of the configuration !•••••) under repeated application of the backward transfer matrix from 
the configuration at i = 0: • oooo) + -- - + |oooo»). 



C. Choice of random maps 



Transition probabilities of the forward transfer matrix must satisfy the Markov chain transition probabilities for 
single particles, but the choice of random maps is otherwise unrestricted. The one-particle sector is trivially correct 
for independent moves as in our diffusion model of Section II A. We now discuss several alternative random maps 
for the one-dimensional diffusion, which may lower (or increase) the coupling time (with, however T coup > r corr ) or 
achieve a rapid reduction of the number of configurations for smaller time scales. 

A naive example for the one-dimensional diffusion example consists of arrows, such as in Fig. 2, but which for one 
time t all point into the same direction, straight, up, and down, each with probability 1/3, so that single-particle 
moves satisfy the detailed balance condition. Evidently, this random map does not couple, and the non-coupling 
Markov chains, in Eq. (8), converge in a time t* = T cor r. We now modify this rigid algorithm by allowing arrows to 
change direction with probability e. This makes the Markov chain couple on a time scale ~ log -, much larger than 
the correlation time r corr , for small e. The choice of independent random moves (e = 1) is optimal in this class of 
maps, but it is not the choice minimizing r coup among all random maps. For example we may choose correlated moves 
for selected neighboring pairs of sites say, for sites (1,2) and (3,4) and let the move from site 5 be independent (see 
Fig. 5). 



Elements of 7 ltw > corr are, for example, 



T fw,corr|| oo .. o )^| oo .. o )} =0 
T fw ' corr {| OO..o) -> | OO.OO)} = 1/3 
T fw ' corr {| OOMO) -> | O • O O.)} = 1/3 
T fw,corr|| oo . #o )^| oo#o# )} = . 
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FIG. 5: Left: Transition probabilities for pairs (1,2) and (3,4) for the correlated random map. Right: In the extended Monte 
Carlo simulation shown, the one-particle transition probabilities are as in Fig. 1, but nearest-neighbor coupling is favored. 

The single-particle sector of this algorithm is as before, but the second-largest eigenvalue of the transfer matrix T fw 
with such correlated pair moves becomes smaller, indicating faster coupling: 

A*' 1 < A' w ' corr = 0.777 < A|, w ' indop = 0.838. 

Coupling times of both "independent-arrow" random mapping and "correlated-pair" random mapping scale alike for 
large N. We note that in applications, as in our patch algorithm of Section V, it might be not so much of interest to 
speed up the coupling than to rapidly decrease the number of possible configurations at times t < T coup . Therefore, 
one goal could be to decrease the eigenvalues of the matrix T k ' k ,k ^> 1, whose time scales correspond to the rapid 
reduction of the number of configurations towards more manageable numbers. 
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D. Exact sampling, coupling from the past 

As discussed in Section II B, the coupling of Markov chains allows one to produce exact samples of the equilibrium 
distribution: In the diffusion example, we were able to run the Monte Carlo simulation backwards in time using T bw , 
but this matrix can usually not be constructed. To find the sample at time t — 0, one may tentatively set a time to < t 
and produce all the random maps between time to and t. One can then check explicitly whether all the possible initial 
conditions at time to have coupled, that is, for the diffusion problem, whether the initial TV-particle configuration 
• ••••} has yielded one of the one-particle configurations. If this goal has not been reached, one must complement 
the random maps already computed with random maps for earlier times (see Fig. 4). The one-dimensional diffusion 
problem without periodic boundary conditions illustrates an algorithm which determines the coupling time with much 
less effort. We consider odd times at which only sites 1, 3 and 5 and even time steps at which only sites 2 and 4 may 
flip (see Fig. 6). This preserves the correct stationary probability distribution, but the trajectories no longer cross 
each other (as at time to + 1 in Fig. 2). As a consequence, it suffices to follow the two extremal configurations, which 
start at sites k = 1 and k — N, from time to on in order to determine the coupling time for a given full Monte Carlo 
simulation. The multiple-particle Monte Carlo simulation starts with the state |»ooo») until it yields a single-particle 
state. The above strategy of following extremal configurations can be applied to the ferromagnetic Ising model (but 
not to spin glasses, see [2, 5]). In this case, the two configurations with all spins up and all spins down, respectively, 
are extremal. This idea also holds for the heat-bath algorithm of two-dimensional directed polymers in a random 
medium. 
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FIG. 6: Extended Monte Carlo simulation with odd (o) and even (e) time steps on a lattice without periodic boundary 
conditions. Trajectories cannot cross, and the coupling of the two extremal initial configurations (the simulations starting at 
time to from sites 1 and 5) determines the coupling time. 
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III. COUPLING AND CONVERGENCE IN SPIN MODELS 

We study the Edwards- Anderson ± J Ising spin glass on a d-dimensional square lattice, where each site is randomly 
coupled to all its 2d neighbors. The energy of a configuration er = (01, . . . , a N ) with <r k — ±1 is: 

{id) 

where indicates the sum over nearest neighbors. This model has a phase transition at finite temperature in 3d 
and at zero temperature in 2d. Sampling spin-glass configurations with Markov-chain algorithms is extremely difficult 
in d = 3 dimensions below the critical temperature, but it is also non-trivial in the two-dimensional case [13]. We 
concentrate here on the study of a local heat-bath Monte Carlo algorithm, for which we apply the coupling-from- 
the-past protocol and obtain exact samples. We note that in two dimensions, the exact partition function of the 
Ising model on a finite lattice can be determined exactly for any choice of couplings [2, 14]. This makes possible a 
direct-sampling algorithm, which is completely unrelated to the material presented here, but which we sketch, for the 
sake of completeness, in Appendix B. 

The heat-bath Monte Carlo algorithm for spin models updates at each time step a randomly chosen site i of a spin 
configuration by comparing a function of the local field on site i with a uniform random number Tj(i) = ran[0, 1]: 

1—1 else 

where hi(t) — J2j Jij^jit) 1S tne local field on site i. A realization of the Markov chain corresponds to sampling the 
real-valued random numbers {. . . , T(i ), T(to+l), . . . , T(— 1)} and the random integers {. . . , i{to), i(t a + l), . . . , i(— 1)}. 
The unit of "physical" time (one "sweep") corresponds to N individual updates. The situation is now much more 
complicated than for the Id diffusion, as the role of the five initial configurations in Fig. 2 is taken up by the 2 N 
possible spin configurations. To prove coupling one must show to which configuration they all converge at the coupling 
time. The state space is huge and one must find strategies to avoid enumerating and surveying 2 N configurations. 

A. Partial-survey approximation 

In [6] , we presented an exact-sampling algorithm which works down to quite low temperatures in the two-dimensional 
Ising spin glass, and which is also operational in three dimensions. We found that practically the same results could be 
obtained by starting the simulation at time to not from all the 2 N initial configurations, but from a more manageable 
number Af(t ) of randomly chosen configurations. We show in Fig. 7 that such "partial survey" calculations yield 
useful lower bounds for the coupling time scale r coup . Each curve in the figure represents the mean number of distinct 
configurations remaining after coupled Monte Carlo simulations (that is, with the same random numbers (T,i) for 
all configurations) for different values of Af(to)- Increasing Af(to) within this partial-survey approximation naturally 
improves the lower bound on the coupling time but, in practice, the value obtained saturates quite quickly. 

B. Correlation functions 

As discussed previously, r coup is always larger than r corr because only non-coupling chains contribute to correlation 
functions (see Eq. (8)). To again illustrate the relation between coupling and convergence times, we separate in 
Fig. 8 non-coupling chains from the calculation of the spin-spin correlation function of a 8 x 8 spin glass at inverse 
temperature (3 = 1. Indeed, even if the chain has not coupled, the configurations at may lose the dependence on the 
initial configuration <r . 

IV. COUPLING FOR HARD-SPHERE SYSTEMS 

In this section we discuss the application of the "Coupling from the past" protocol to hard-sphere systems. The 
study of Monte Carlo algorithms for hard-sphere systems goes back a long time, as the Metropolis algorithm was 
first implemented for hard disks, that is, two-dimensional spheres [1]. Even today, the physics of the hard-disk system 
is not well understood, and Monte Carlo algorithms have not been developed as successfully as, say, for the Ising 
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FIG. 7: Number of configurations of the partial survey approximation with J\f(to) random initial configurations for the Ising 
spin glass on a 16 x 16 lattice at temperature /3 = 0.5. We average over 10 choices of Jij, and use the same values of Jij and 
the same random numbers T, , i for all initial configurations. 
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FIG. 8: Spin-spin correlation function for all Markov chains and for the non-coupling Markov chains only (/? 
initial conditions) (Two-dimensional 8x8 spin glass at f3 — 1). 
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model. In this very constrained system, the estimation of correlation times is quite controversial, especially at high 
densities [18], and rigorous results from exact-sampling approaches would be extremely welcome. 

We first discuss the birth-death formulation of the Markov-chain Monte Carlo algorithm for this system and then 
compute lower bounds on coupling times using the partial-survey algorithm. Its empirical coupling time saturates 
(for increasing jV(to)) to much smaller values than the coupling times obtained by the summary-state method [9, 10]. 
This suggests that these previous algorithms are not optimal, an impression which is confirmed by our local-patch 
algorithm of Section VE. 

The partition function of hard spheres in the grand-canonical ensemble, with fugacity A, is given by a weighted sum 
over legal configurations of spheres: 

+oo „ 

2 = E d 2N ° (N) A" 6 (aW) . (11) 

Here, configurations of N spheres are written as: 

o-W = {(xi,yi), (x 2 ,y 2 ), (x N ,y N )}, 

where (xk, yk) denotes the centers of the spheres. In Eq. (11), Q(a^ N ^) equals one if spheres of the configuration (jW 
do no overlap and zero otherwise. We again use periodic boundary conditions. 
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A. Birth— death algorithm for hard spheres 

The spatial Poisson birth-death process allows us to apply coupling from the past to hard-sphere systems (see 
[9, 10]): Disks of radius r arrive ("are born") randomly on a two-dimensional unit square with constant rate A. Once 
born, they disappear ("die") with unit rate. 

The probability for a disk to arrive within an infinitesimal time dt in a small box of area dS centered at (x, y) is 
XdtdS. This disk is added to the configuration only if it overlaps with no other disk present. Each disk disappears 
with probability dt within the time interval dt. Sphere are added at point (x, y) or removed from the configuration 
according to the detailed balance condition. With the notation a^ N+1 ' ) = o-W U {(x,y)} we have: 

P(/l - <j (N +^)n(a^) = A Q(<j (N+ ^)n(a^) 

= 7r(a( Ar + 1 )), 

7V (JV+1) - ff (S) )^ (f,+1) ) = 1 x 7T(a( JV + 1 )). 

In Fig. 9, we illustrate the time-evolution of accepted and rejected birth-death events on a one-dimensional hard- 
sphere problem, starting from an empty initial condition at time to- 




FIG. 9: Simulation of the birth-death algorithm for one-dimensional "spheres" in a box. The simulation starts at time to and 
stops at to + tsim with N — 2 spheres. Transparent spheres are rejected because they overlap with spheres already present. 

In the hard-sphere algorithm, the probability distribution of time intervals between successive births is an exponen- 
tial with parameter A: V(i~b) = Ae~ Arb . In Fig. 9, the life time of a sphere is represented by a horizontal extension of 
the box, irrespective of whether it has been accepted or not (the vertical dimension denotes the diameter) . Life times 
are exponentially distributed as well. For the exponential distribution, the time before the next death of a system 
of N spheres follows an exponential distribution with parameter N. Likewise the time before any event, (birth or 
death), follows an exponential distribution with parameter A + N. The probability for the next event to be a birth is 

thcn XTF- 

B. Coupling and partial survey approximation 

Coupling from the past applies to hard-sphere systems even though the space of configurations is continuous (unlike 
in lattice simulations). To apply the protocol, one considers a time evolution, as in Fig. 9, but stretching back to time 
t = — oo. Two special aspects must now be handled: 

First, we must determine which boxes (corresponding to spheres) are indeed placed ("True"), and which ones are 
rejected ("False"). This is difficult to decide at high density. However, in the low-density case presented in Fig. 10, 
several spheres are "True" , simply because they do not overlap with already present "True" or "False" spheres. This 
allows the status of other spheres to be fixed and, finally, the configuration to be constructed. In the limit N — ► oo, 
the approach works up to a constant density [17]. This density is much higher than the density oc 1/N direct-sampling 
algorithm can achieve [2]. This approach [9, 10] is equivalent to deciding whether a given spin is up or down in the 
"summary state" algorithm for Ising systems [7, 8], which, in the thermodynamic limit works down to a fixed constant 
temperature. 

Second, one must fix the initial condition at time —T, because spheres born at times smaller than — T may still be 
alive at time — T . This is solved through the sampling of a second time, T star t, after which we know that all spheres 
present at time — T have disappeared. The time interval T sta rt + T is sampled as the maximum of N max life times, 
where 7V max is an upper bound on the number of spheres in a legal configuration. 

Figure 10 sketches the time evolution of a Monte Carlo simulation for the one-dimensional hard-sphere problem 
which has started at time t — -co. Boxes are drawn starting from time — T, but the simulation is picked up at time 
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FIG. 10: Time evolution of a one-dimensional birth-death simulation in a box of size L. All spheres correspond to rectangles 
whose horizontal extension indicates their life time. From any possible cut of "True" boxes at r sta rt (7 boxes actually cut the 
line, so there are < 2 7 possible cuts) one can deduce the output at time t — 0, as in Fig. 9. 

Tstart- It is straightforward to complement the simulation shown (between times —IT and — T, for example), in case 
it does not couple in the interval shown. However, we must show that it couples between — 2T and — T or at least 
results in less than iV max spheres. In the Monte Carlo simulation in Fig. 10, the status of the boxes at later times can 
be easily decided, because at later times all spheres belong to clusters which are disjoint from the initial condition. 
However, this possibility disappears at higher densities. A simple example of this is shown in Fig. 11. As one cannot 
decide on the status of the initial sphere (which crosses the line at T start ), we should initialize the simulation with the 
two configurations, one corresponding to a "True" state and one to a "False" state. After several steps of the time 
evolution, we arrive in both cases at the same physical configuration (the two dark spheres, which are both "True"). 




— i — 

start 



FIG. 11: Example of a time evolution of the one-dimensional birth-death simulation which where no single sphere can be 
decided independently. Starting with all possible choices of the initial configuration at t = T sta rt allows to prove coupling. (The 
two dark spheres are "True", while the transparent sphere must be "False".) 

For all times t > T sta rt; we consider the set C{t) of all "True" or "False" spheres crossing the time line at t (see 
Fig. 10). From the set C, one can in principle construct all the possible initial configurations, but their number remains 
huge. As in the spin-glass case, we may also select Af(T stalt ) among these configurations, and propagate these. This is 
again the partial survey approximation. In Fig. 12 we compare average lower bounds on the coupling time from this 
approximation with results from the summary state algorithm[9]. In the time-evolution of Fig. 10, one can determine 
the number of remaining configurations at any time t > T star t and detect when exactly the coupling occurs. 

V. LOCAL-PATCH ALGORITHM 

In the present section, we discuss our local-patch algorithm, which performs the heat-bath dynamics for a general 
(i-dimensional Ising spin glass on an A^-site hyper-cubic lattice. This algorithm allows us to control all the 2 N initial 
configurations even for very large lattices and to eventually prove that the system has coupled. The Python script 
implementing this algorithm has less than 300 lines. It is available electronically and a listing of the code is contained 
in Appendix C. 

A. Patches 

The (non-rigorous) partial-survey algorithm of Section III A determines the coupling time for a subset of all the 
configurations at time t. The (rigorous) patch algorithm, in contrast, works with a superset of all configurations at 
time t: by restricting the configurations to the smaller region of a patch, one severely limits their number, at the price 
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FIG. 12: Coupling times of the summary-state algorithm for hard disks [9, 10] compared to the local-patch algorithm. Lower 
bounds are provided by the partial-survey approximation. The disks' radius is r — 0.04 in a unit square box with periodic 
boundary conditions, so that there are ~ 60 disks at density rj = 0.3. 



of introducing compatibility problems between neighboring patches. For the two-dimensional spin glass, we use N 
rectangular patches of same shape and orientation, with M sites, and initially at t — to, we have 2 M spin configurations 
on each patch. Likewise, the set of global spin configurations is broken up into a list [Si(t), . . . , S'jv(t)] of sets Sk(t) 
of spin configurations restricted to patches k. We can recover a superset of all relevant spin configurations from 
the direct product 

n(t) = Si(t) <8> S 2 (t) ® ■ • ■ ® Sjv(i)/(compat). (12) 

Here, each configuration of fi(t) is pieced together from configurations on all patches, with compatible spins on all 
lattice sites. (Two compatible spin configurations, on patches k and I, are shown in Fig. 13). On large lattices, the 
direct product in Eq. (12) can be performed only if the number of spin configurations per patch is small. If there is 
only one configuration per patch, we can construct a unique global configuration on the whole lattice. 

For each time step t of the heat-bath algorithm, we choose a random lattice site i and a random number T = ran[0, 1] 
and then update the spin at for all configurations on all patches containing i. The site i may be in the center of a 
patch k (all the neighbors of site i also belong to k, as in Fig. 13). In this case, each configuration of Sk(t) yields one 
configuration of Sk(t + 1). Several configurations in Sfc(t) may yield the same configuration in Sk{t + 1), so that the 
size of Sk does not increase in this case. If the site i is on the boundary of a patch /, we only know upper and lower 
bounds for the field on the site i and, depending on the value of the random number T, may be unable to update Oi. 
In this case, we add two configurations to the set Skit + 1)> corresponding to at = —1 as well as <7i = +1. The set 
Sk(t + 1) may then contain more configurations than Sfc(i). 



+ - 

+ - + + 
+ + + - 



patch k spin configs patch I 

FIG. 13: Two patches, k and I, with a pair of compatible spin configurations. 
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B. Compatibilities, pruning 



Besides updating configurations on patches, we also perform a "pruning" operation: Figure 13 presents two "com- 
patible" configurations on patches k and I. These could possibly be pieced together into a global configuration, 
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together with configurations on other patches. On the other hand, if the set <S) contains no configuration compatible 
with a configuration er on patch k, we can eliminate (prune) <x from S^. Pruning may be implemented through a 
dictionary (hash table), using as "key" the part of the patch configuration in the overlap region between k and I, 
and as "value" the list of patch configurations sharing this key (see Fig. 14). This is programmed very easily in the 
Python programming language (see Appendix C). 



key 



value 

[EH, 













* 














El 




* 












:;; 










5 


*** 












> 







Hz: 



FIG. 14: A dictionary (hash table) associating keys (configurations in the overlap region between patches k and I) to values 
(lists of patch configurations with the given key in patch k). 



Pruning can be iterated until all the sets Sk(t) are pairwise compatible. To achieve this goal, it suffices to prune 
nearest-neighbor patches only. We have found it useful to perform one pruning operation for each pair of nearest- 
neighbor patches after a certain number of updates (see Fig. 15). The average number of configurations per patch 
saturates to a value which depends on the temperature and also the size of the patch. This is due to the balance 
between the decrease of the number of configurations induced by the coupling and its increase caused by the noise at 
the patch boundaries. This noise is reduced through the crucial pruning step of the algorithm. 




FIG. 15: Average number of configurations per patch vs. time, at inverse temperature /3 = 1 in the two-dimensional Ising spin 
glass on a 16 x 16 lattice at temperature (3 = 1. For a given patch size (here 3x3) the number of configurations per patch 
saturates to a value that depends on the number of prunings per sweep (here: from one to eight). Each pruning is done once 
for all nearest-neighbor patches. 



C. Merging of patches 



As illustrated in Fig. 15, the number of configurations per patch does not necessarily drop to one at large times, 
even if the underlying heat-bath dynamics couples. The entropy per spin is smaller for larger patches, because the 
influence of the boundaries is reduced. However, one cannot start the computation with large patches because of the 
large number of possible configurations. A merging procedure allows to increase the patch size in a rigorous way. 
Merging is implemented analogously to pruning: for overlapping patches k and I, dictionaries are again computed 
with the same keys and values (see again Fig. 13). For a given key configuration on the overlap region, we assemble 
the corresponding values in the dictionary of patch k with all corresponding values in the dictionary of patch I. All 
these couples of configuration must be taken into account for the larger patch kUl. The merging of the configurations 
on neighboring patches can be implemented very efficiently in the Python programming language (see Appendix C). 
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In our computations, we start with small square patches, say of size 3x3, and then pass to the size 4x3, after a few 
sweeps, then to size 4x4, etc. An analogous procedure is followed in higher dimensions. Results obtained with this 
"jump-start" approach are shown in Fig. 16 for the two-dimensional Ising spin glass at temperature (3 — 0.5, with a 
disorder average performed over about 100 samples. 




sweeps 

FIG. 16: Number of configurations per patch for the 32 x 32 spin glass at temperature (3 = 0.5. A simulation with constant 
3x3 patches is compared to the result of a "jump-start" procedure 3x3^4x3^4x4->---^5x5. All results are 
averaged over 100 samples. 



D. Memory of compatibilities, variants 

In the patch algorithm, the pruning procedure detects inconsistent configurations in a particular patch. Two 
configurations on different patches are considered compatible if their spins match in the overlap region (at time t). 
More generally, we can keep track of the past evolution of patch configurations and may then declare them compatible 
only if they have matched for all times up to t. Otherwise, they cannot belong to a unique global spin configuration. 

In [6], the bipartite nature of the square lattice was used to update one entire sub-lattice at a time. In this approach, 
only one sub-lattice is stored at a time. This allows one to start with larger patches, but the compatibilities between 
configurations are less well conserved. By contrast, in the present algorithm, we keep the information on both sub- 
lattices, and one of the sub-lattices is the past configuration. Two configurations are compatible in this new version 
if they are compatible on both sub-lattices. Likewise one can use past values of a configuration a(t) to restrict its 
compatibilities with configurations on other patches. 

Other generalizations are more straightforward, one can for example optimize the shape of patches in order to 
minimize the number of spins on the boundary, and work with more than N patches in order to increase the chance 
for detecting incompatibilities of configurations. 

E. Exact sampling for hard spheres with local patches 




patch k disk configs patch / 



FIG. 17: Breaking up disk configurations into patches, with two patches k and I shown. 
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In this section, we adapt the local-patch algorithm for the classical model of hard disks in a two-dimensional box 
with periodic boundary conditions. As mentioned, it works for large system sizes at higher density than previous 
method of exact sampling. 

In Section IV B, we introduced the set C(t) of all disks the dynamics has tried to add in the box and which have not 
disappeared at time t. The coupled Monte Carlo simulations start with all possible configurations that are allowed 
by the set C at time T sta rt- Because of possible overlap between disks some of the 2 #c possible configurations arc 
invalid but there may be far too many of them in practice. To reduce the number of configurations which must be 
handled, we introduce a regular square lattice with N sites covering the simulation box. Likewise, a superset of all 
feasible configurations on a patch at time T sta rt can be deduced from the set C, restricted to disks (True or False) 
with centers inside the patch. From then on, whenever disks appear in the simulation box, we can decide whether 
they are accepted on a particular patch configuration by checking overlaps into the patch only. Disks that disappear 
are simply removed from all the concerned patch configurations. At birth time, if the disk to be placed on a patch 
configuration may overlap with a disk outside the patch, the configuration is split into two: one configuration with 
the new disk and one without (as for spin systems). To detect and prune irrelevant configurations we check that 
the updated sets of configurations are compatible with other patches by a pruning procedure analogous to the one of 
Section VB. After several updates, the pruning is performed for most of overlapping patches several times. 

key value in patch k 

s -[<ta.4g].<&<©] 

key value in patch 1 

B Mi 

FIG. 18: The pruning of a pair of patches via the construction of dictionaries. Each dictionary associates keys (configurations 
in the overlap region) to values (lists of patch configurations with the given key) (compare with Fig. 14). The merging of 
patches k and I would lead in this example to 4 x 3 = 12 configurations on the combined patch. 

For any pair of overlapping patches k and /, the pruning eliminates patch configurations with are inconsistent with 
all other configurations on a neighboring patch. As in the spin glass case, this process can be implemented with 
dictionaries (hash tables) (see Fig. 18) and can be used to merge configurations on neighboring patches into larger 
local configurations. 

Figure 12 displays the mean coupling times of a hard-disk birth-death simulation for several choices of the fugacity 
A. The radius of the disks is r = 0.04 in a unit square box with periodic boundary conditions. Results are compared 
to the partial-survey approximation algorithm, with A/"(T start ) = 1000 initial configurations and to the results of the 
summary-state algorithm. 

We concentrated on determining the coupling times of the birth-death dynamics for hard disks. However, the 
regime of operation of this algorithm is far in the liquid phase (see Fig. 12), and the physically interesting regime, 
around the liquid-solid transition density, r\ ~ 0.71 [2] is still out of reach for exact sampling methods. For hard disks, 
it remains a challenge to set up a working partial-survey algorithm with correlation times comparable to those of the 
usual Metropolis algorithm [18]. 

VI. CONCLUSION 

In this paper, we have discussed exact-sampling algorithms which allow one to totally eliminate the influence of 
the initial condition from a Markov-chain Monte Carlo simulation. This overcomes one of the main limitations of the 
method, namely the rigorous estimation of the correlation time. We discussed central subjects, such as the relation 
between coupling times and convergence times, in a simple example of one-dimensional diffusion, before applying 
them to Ising spin glasses and to hard-sphere simulations. Algorithmically, the exact-sampling framework obliges 
one to follow the entire state space of a system. In the absence of simplifications, such as the half-order discussed in 
Section II, this can be done approximately through a partial survey of A/"(*o) initial conditions. One can also restrict 
the configuration in size onto so-called patches, thereby restricting their number. A superset of the set of global 
configurations can in principle be reconstituted from the patches. This is easier when the patches are large, and we 
showed how pruning and merging operations allow one to increase the size of patches during the simulation and to 
finally prove coupling. Our exact-sampling algorithm works both for spin glasses and hard-disk systems, and were 
able to go to lower temperatures, and higher densities than previous methods. 
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The partial survey algorithm, which can be implemented easily, allowed us to prove that our local-patch algorithm 
is optimal for the local dynamics for both spin-glass and hard-disk systems. We have provided a number of new idea 
in order to allow exact-sampling methods to reach the phases transitions of the three-dimensional Ising spin glass and 
the critical density of hard disks. 

APPENDIX A: SIMILARITY OF FORWARD AND BACKWARD TRANSFER MATRICES 

Forward and backward transfer matrices completely describe the coupling dynamic of general Markov chains (on 
a finite state space), that is their elements are the coupled transition probabilities between sets of configurations. 
Forward and backward matrices represent "extended" Monte Carlo dynamics, in the two time directions. These two 
formulations are equivalent, even though the forward dynamics, starting from t — 0, does not generate exact sam- 
ples. Here, we demonstrate similarity between forward and backward transfer matrices, and construct the similarity 
transformation between the two. 

Let fi be the finite space of configurations of the problem of interest, may contain all TV positions on the iV-site 
diffusion problem or the 2 N configurations of a TV-site spin system. 

The /c-configuration states build up an "extended" state space. They provide a natural basis for the forward and 
backward matrices. This basis is 2° — 0, the set of non trivial parts of 0. For any state I e 2° — we define / as the 
set of states of the basis J that has at least one configuration in common with I (Je/o/nJ^8). The similarity 
matrix P is then defined as: 

P\I):=-£\J) (Al) 

Jei 

A random mapping — arrows for the case of ID-diffusion — is a mapping on £1: / : O — » Q. It defines a time step of 
the Markov chain for every configuration and satisfies V(f(i) = j) = p(i — > j) and its weight is noted w(f). In the 
case of Id-diffusion with independent arrows, or any "independent" random map in general, we naturally define the 
weight of the random map as a product of elements of the Monte Carlo transfer matrix as 

w(f) = U ienP (i - /(»)). 
The forward matrix associates any state I to all states that it is connected to by a mapping: 

T tw |7) = w (fW))- 

f rand, map 

Using Eq. (Al) we find 

PT fw |/)= £ w{f) £ |J). (A2) 

/ rand, map J£f(I) 

The backward matrix T bw has different rules but we will show that the similarity P T iw = T hw P holds. Using a 
random map /, a state J at time t evolves to another state K at time t + 1 in the backward process if and only if 
/ _1 (.ff) = J. For example in the ID-diffusion a hole goes to a hole and a particle goes to a particle. Therefore: 

/ rand, map K,f- 1 (K) = J 

and finally: 

T bw P|/)= £ h/)E E \ K )= E E i*>- ( A3 ) 

/ rand, map jgj K,f~ 1 (K)=J f rand, map KJ- 1 (K)eI 

In fact Eqs (A3) and (A2) are equivalent because / -1 (.?r) overlaps I if and only if K overlaps /(/) (/ _1 (if) P\ I ^ 
K (1 f(I) 7^ 0). This proves the similarity of the backward and forward matrices. 
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APPENDIX B: EXACT SAMPLING OF TWO-DIMENSIONAL SPIN GLASS USING ANALYTIC 

SOLUTION. 



In this appendix, we sketch for completeness an unrelated direct-sampling algorithm for the two-dimensional Ising 
spin glass. To generate exact samples, this algorithm does not use Markov chains. It rather relies on the fact that the 
partition function of the two-dimensional Ising model or of one sample of the spin glass on a planar lattice with N sites 
can be expressed as the square root of the determinant of one AN x AN matrix (for open boundary conditions) or of 
four such matrices (for periodic boundary conditions) [14]. This relation has been much used in the recent literature, 
in order to study the physics of the two-dimensional Ising spin glass at low temperature [15, 16]. The partition 
function yields the thermodynamics of the system, but the knowledge of entire configurations gives for example access 
to complicated spatial configuration functions. 

The sampling algorithm for two-dimensional spin-glass configurations constructs the sample one site after another. 
Let us suppose that the gray spins in the left panel of Fig. 19 are already fixed, as shown. We can now set a fictitious 
coupling J* ; either to — oo or two +oo and recalculate the partition function Z± with both choices. The statistical 
weight of all configurations in the original partition function with spin "+" is then given by 

= Z+cxp (pj kl ) 

+ Z+exp (/?J fcJ ) + Z_exp(-/3J fcJ )" 1 ; 

and this two-valued distribution can be sampled with one random number. Equation (Bl) resembles the heat-bath 
algorithm of Eq. (10), but it is not part of a Markov chain: After obtaining the value of the spin on site k, we keep 
the fictitious coupling, and add more sites. Going over all sites, we can generate direct spin-glass samples at any 
temperature. We note that this algorithm is polynomial, and the effort is basically temperature-independent, both 
for the two-dimensional Ising model and the Ising spin glass (see also [19]). 
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FIG. 19: One iteration in the direct-sampling algorithm for the two-dimensional Ising model. The probabilities ixiou = +1) 
and 7r(<7 fe = —1) (with k the central spin) are obtained from the exact solution of the Ising model with fictitious couplings 
Jjk = ±oo. 



APPENDIX C: LISTING OF PYTHON CODE 



The following Python code has produced all the spin-glass data presented in this article. An analogous program 
was used for the hard-disk system. An electronic version of the code is available from the authors. 

Algorithm 1: pruning-ND.py 



#!/ usr / bin /python 
## 



## PROGRAM 

## PURPOSE 

## 

## 

## 

## OUTPUT 
## 

## VERSION 
## AUTHOR 
## LANGUAGE. 
## 



pruning-ND . py 

This program performs the heat— bath and the pruning on a 
N— dimensional hypercubic lattice with periodic boundary 
conditions and with cuboid patches. 
Version with jump— start capability . 
mean number of configurations per patch vs. time 
(can be modified) 
04-OCT-2009 
W. Krauth , C. Chanal 
Python 2.5 



from random import uniform , randint , seed 
from operator import itemgetter 



shuffle , choice 



16 
17 
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19 
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23 
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32 
33 
34 
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39 
40 
41 
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import time, math, os , sys 
## 



## Sample geometry 
##- 



def t or us.ne igh bor s ( N_dim , L ) : 
N = L**N_dim 
site_dic = {} 
coord_dic = {} 
for j in range(N): 
j-d = j 

x.d = j//(L**(N_dim-l)) 
coord = [x_d] 

for d in range ( N_dim — 1 ) : # this loop does not run anything for N_dim 
j_d = ( j_d— x_d*L**(N_dim— d — 1)) 
x.d = j_d//(L**(N_dim-d-2)) 
coord . append ( x_d ) 

coord . reverse ( ) # optional , to set the usual order of the directions 
coord=tuple (coord) 
s i t c _d i c [ coord ] = j 
coord_dic[j] = coord 
nbr = [] 

for j in range (N): 

coord = 1 i s t ( coor d_di c [ j ] ) 
nbr_list = [ 
for d in range ( N_dim ) : 

coord [d] = (coord [d] + l)%L 
coord_p = tuple(coord) 
n b r _1 i s t . append (sitc_dic [coord -P ] ) 
coord [d] = (coord [d]-l+L)%L 
for d in range ( N_dim ) : 

coord [d] = (coord [d]-l+L)%L 
coord.p = tuple(coord) 
nbr _1 i s t . append (sitc.dic [coord.p]) 
coord [d] = (coord [d] + l)%L 
nbr.list = tuple (nbr.list) 
nbr . append ( nbr.list ) 
nbr = tuple ( nbr ) 
return nbr , site.dic , coord.dic 
## 



## Initial patch geometry 
##- 



def pat ch_set ( N_dim , L ,M, s i t e _d i c ,coord_dic): 
patch = [j 
for j in range (N): 

coord = 1 i s t ( coor d_di c [ j ] ) 
dummy_list = [ 
for k in range (M**N_dim ) : 
coord_p = [ 

for d in range (N_dim ) : 
1 = (k//(M**d))%M 
coord.p . append (( coord [d]+ 1)%L) 
coord.p = tuple ( coord.p ) 
dummy.list . append (sitc.dic [coord.p]) 
patch . append ( tuple ( dummy.list ) ) 
patch = tuple (patch) 
return patch 
## 



## Neighbor relations on patches 
##- 



def nbr_patch_set( example .patch, nbr): 
nbr_patch = [ 
for j in example_patch : 
dummy = [ 
for k in nbr [ j ] : 

if k in cxample_patch : dummy . append ( example_patch . index ( k ) ) 
dummy = tuple (dummy) 
nbr_patch . append (dummy) 
nbr_patch = tuple ( nbr_patch ) 
return nbr.patch 

## 
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## Initial patch configurations 
## : 



def p a t c h _i n i t ( patch ) : 

def bin (n , conf .length ) : 

## 

## convert n to binary number with conf .length digits 
## 

q = -1 
bin_conf = 
n.digits = 
while q != 0: 

q = n//2 
r = n%2 

bin_conf = 'r'+bin_conf 
n = q 

n.digits = n_digits+l 

n.digits = conf.length — n.digits 

for i in range ( n _d igit s ) : 
bin_conf = '0' + bin_conf 

return bin_conf 
configs = [ 
dummy_list = [ 
number = len ( patch [ ] ) 
for j in range (2 ** number ) : 

x = bin(j , number) 

dummy_list . append (x) 
for i in range ( len ( patch )) : 

configs . append (set (dummy _list ) ) 
return configs 



## 

## Updating configurations on all patches 

##- 



def confs.update (N_dim , configs, nbr_p at ch,Jij ,patch,i_site ,Upsilon): 
for k in range (N): 

if i_site in patch [k]: 
dummy = list(patch[k]) 
pos=dummy . indcx( i_sitc ) 

Hist = [patch[k][x] for x in nbr.patch [ pos ] ] 
Jij.list = [Jij [ ( i _si t e ,m)] for m in Hist] 
config_k = configs [k] 
config_kp = set ([]) 
for c in config_k : 

ficld_char = itemgetter (* nbr_patch [ pos ] ) ( c ) 

ficld_sum = sum ( ( 2 * eval (m) — 1)* J for (m,J) in zip ( f icld_char , Jij.list )) 
field_min = field_sum — 2*N_dim+len ( f ield_char ) 
field_max = field_sum+2*N_dim— len ( f ield_char ) 

b.onc = c [ : pos] + ' 1 '+c [ pos + 1 : ] # co nfiguration with ' 1 ' at position 'pos ' 
b_zcro = c [ : pos]+ ' '+c [ pos + 1 :] # configuration with '0' at position 'pos 
if Upsilon < phplus [ field_min ] : 

config_kp .add(b_one) 
elif Upsilon > phplus [ field_max ] : 

config.kp . add (b.zcro) 
else : 

config.kp . add (b.zcro ) 
config.kp . add ( b.one ) 
configs [k] = config.kp 
return configs 

## : 



## Pruning of two patches 
##- 



def pr une_pair ( conf ig _k , co nf i g _1 , KEY_k , KEY_1 ) : 
f = itemgetter (*KEY_k) 
Dick = {} 
for x in config_k : 

a = f(x) 

if Dic_k . has_key ( a ) : Dic_k [ a ] . add (x ) 
else: Dic_k[a] = set ( [x] ) 
sct_k = set ( Dic_k . keys () ) 
f = itemgetter (*KEY_1) 
Dic.l = {} 



for x in config_l: 
a = f(x) 

if D i c _1 . has_key ( a ) : Dic_l [ a ] . add (x) 

else: Dic_l[a] = set ( [x] ) 
set_l = set ( Dic_l . keys () ) 
set_kl = set.interscction(set_k,sct_l) 
conf ig _k = set ( ) 
config_l = set() 
for x in se t _kl : 

config_k . update (Dic_k [x ] ) 

config_l .updatc(Dic_l[x]) 
return config_k , config_l 



## 

## Merging of two patches 

##- 



def merge.pair ( config.k , config.l , KEY_k , KEY J , KEY.add ) : 
f = itemgetter (*KEY_k) 
Dic.k = {} 
for x in config.k : 

a = f(x) 

if Dic_k . has_kcy ( a ) : Dic_k [ a ] . add (x ) 

else: Dic_k[a] = set ( [x] ) 
sct_k = set ( Dic_k . keys () ) 
f = itemgetter (*KEY_1) 
Dic.l = {} 
for x in config.l: 

a = f(x) 

if D i c _1 . has_key ( a ) : D ic _1 [ a ] . add (x ) 
else: Dic_l[a] = set ( [x] ) 
set_l = set ( Dic_l . keys () ) 
set_kl = sct.intcrscction(sct_k,set_l) 
config.kl = set () #merged set 
for x in se t _kl : 

for y in Dic_k[x]: 

for z in D ic _1 [ x ] : 

zprime = ' ' . j oi n ( i t emget t er ( *KEY_add ) ( z ) ) 
config.kl . add (y+zprime) 
return config_kl 
## 



## main program starts here 



## — 

#seec 
beta 
## 


1(13) 

=0.25 # inverse of 


temperature 




## h 
## 


iat bath definitions 


(see SMAC Fig. 5.20, and SMAC eq . 


5.18) 



N.dim = 3 
phplus = {} 

for d in range ( N_dim + 1 ) : 
field = 2*d 

phplus [ fi e 1 d ] = 1/ (1 + math . exp ( — 2* f i e 1 d * beta ) ) 
phplus [— field ] = 1/ (1 + math . exp (2 * f i e 1 d * beta ) ) 

## 

## loop over samples starts here 

## 



L = 6 

os.systcm( ' echo „ ' hostname ' „ ' date ' ' ) 
print beta, L, '_beta_L' 
for nsamp in range (100): 
N = L**N_dim 

M = 2 #M~N_dim is the initial size of patches 

ovcr.min = (M— 1)*M** ( N_dim — 1) # minimum overlap for the pruning 
dcl_t = 40 # time lap on each patch size 
n.gen = 7 # number of generations 
t_max = 2000 # total number of sweeps 

N_frac = 6 # do N/ N_frac spin updates between prunings 
nbr,site_dic , coord_dic = t or us _nci ghbor s ( N_dim , L) 



## 

## Jij : a dictionary (i,j) — > J-ij 
## ^ 



226 
227 
228 
229 
230 
231 
232 
233 
234 
235 
236 
237 
238 
239 
240 
241 
242 
243 
244 
245 
246 
247 
248 
249 
250 
251 
252 
253 
254 
255 
256 
257 
258 
259 
260 
261 
262 
263 
264 
265 
266 
267 
268 
269 
270 
271 
272 
273 
274 
275 
276 
277 
278 
279 
280 
281 
2S2 
283 
284 
285 
286 
287 
288 
289 
290 
291 
292 
293 
294 
295 



Jij = {} 

for k in range (N): 

for d in range ( N_dim ) : 

Jij [(k, nbr [k] [d] ) ] = choice ([- 1 , 1] ) 

Jij [(nbr[k][d] ,k)] = Jij [(k,nbr[k] [d])] 
patch = p at ch_set ( N_dim , L ,M, s i t c _d i c ,coord_dic) 
nbr_patch = nbr_patch_sct( list (patch [0]), nbr ) 
permut = [k for k in range (N)] 
configs = pat ch _init ( patch ) 

## 

## Loop over generations in the jump— start procedure 

## 

for itcrl in range ( n.gen ) : 
i_dir = itcrl % N_dim 
tot.it = 0. 

patch.size = len ( patch [ ] ) 

print patch.sizc , ' „ s i z e _ of _ pa t c h 

while (tot_it < dcl_t and itcrl < n_gcn — 1) or (itcrl = n.gen— 1 and \ 
tot_it +dcl_t * iter 1 < t.max ) : 
quality = sum ([len(configs[k]) for k in range(N)])/float(N) 
sys . stdout . flush () 

print t o t _i t +d e 1 _t * it er 1 , quality 

#if ( it er 1 ==n_gen — 1 ) : N_Jrac=18 # do N/N_frac spin updates between prunings 
for itcr3 in range (N/ N_frac ) : 

tot.it += l./N 

i_site = randint (0 ,N— 1) 

Upsilon = uniform(0,l) 

configs = confs.update (N_dim , configs, nbr .patch, Jij ,patch,i_site , Upsilon) 

## ■ ' 

## Pruning starts here 

## 

shuffle (permut) 
for kk in range (N): 
k = permut [ kk ] 
for 11 in range (kk + l,N) : 
1 = permut [11] 

inter.set = s et ( patch [ k ] ) & set ( patch [ 1 ] ) 

if len ( inter.set ) >= ovcr.min: # minimum overlap 

KEY_k = [ 1 i s t ( patch [ k ]). index ( i ) for i in inter.set] 
KEY_1 = [ 1 i s t ( patch [ 1 ]). index ( i ) for i in inter.set] 

configs [k] , configs [1] = prune_pair(configs [k] , configs [1] , KEYJc , KEY_1) 
if ( quality = = 1): break 
if ( quality ==l):break 



##- 



## Merging starts here: producing new patches , and computing configurations on them 

##- 



patch_new = [ 
configs_new = [ 
for k in range (N) : 
KEY.add = [] 

dummy _new = list(patch[k]) 

1 = nbr [ k ] [ i _d i r ] # I is the neighbor of k in the 'i-dir' direction 
for x in patch [ 1 ] : 

if x not in patch [k]: 
dummy _ncw . append (x) 

KEY_add . append ( 1 i s t ( patch [ 1 ] ) . index (x ) ) 
patch.new . append (tuple (dummy .new ) ) 
inter.set = s ct ( patch [ k ] ) & set ( patch [ 1 ] ) 
KEYJc = [list (pat ch[k]).index(i) for i in inter.set] 
KEY_1 = [list(patch[l]).indcx(i) for i in inter.set] 

config_k_merge = merge.pair (configs [k], configs [1], KEY_k, KEYJ , KEY_add) 
configs.new . append (config.k .merge) 



## 

## Merging (final steps ): computing neighbor relations on the new patches 

## : 



patch = tuple ( patch_new ) 

nbr .patch = nbr_patch_set(list(patch[0]), nbr ) 
configs = configs.new 
if i _ d i r = N_dim — 1 : 
M += 1 
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296 1 ovcr.min = (M— l)*M**(N_dim — 1) 
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